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Abstract 

Under certain circumstances, Voigt waves can propagate in a biaxial composite 
medium even though the component material phases individually do not support Voigt 
wave propagation. This phenomenon is considered within the context of the strong- 
permittivity-fluctuation theory. A generalized implementation of the theory is devel- 
oped in order to explore the propagation of Voigt waves in any direction. It is shown 
that the correlation length — a parameter characterizing the distributional statistics 
of the component material phases — plays a crucial role in facilitating the propagation 
of Voigt waves in the homogenized composite medium. 
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1 Introduction 

A defining characteristic of metamaterials is that they exhibit behaviour which is not exhib- 
ited by their component phases [1]. A prime example is provided by homogenized composite 
mediums (HCMs) which support Voigt wave propagation despite their component material 
phases not doing so [2]. Although they were discovered over 100 years ago [3], Voigt waves 
are not widely known in the optics/electromagnetics community. However, they have re- 
cently become the subject of renewed interest [4, 5], and the more so in light of advances in 
complex composite mediums [6]. 

A Voigt wave is an anomalous plane wave which can develop in certain anisotropic medi- 
ums when the associated propagation matrix is not diagonalizable [7] . The unusual property 
of a Voigt wave is that its amplitude is linearly dependent upon the propagation distance 
[8]. 

In a recent study, the Maxwell Garnett and Bruggeman homogenization formalisms were 
applied to show that Voigt waves may propagate in biaxial HCMs provided that their 
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component material phases are inherently dissipative [2]. However, the Maxwell Garnett 
and Bruggeman formalisms — like many widely-used homogenization formalisms [9] — do 
not take account of coherent scattering losses. The strong-permittivity-fiuctuation theory 
(SPFT) provides an alternative approach in which a more comprehensive description of the 
distributional statistics of the component material phases is accommodated [10, 11]. In the 
bilocally-approximated implementation of the SPFT, a two-point covariance function and 
its associated correlation length L characterize the component phase distributions. Coherent 
interactions between pairs of scattering centres within a region of linear dimensions L are 
thereby considered in the SPFT, but scattering centres separated by distances much greater 
than L are assumed to act independently. Thus, the SPFT provides an estimation of coher- 
ent scattering losses, unhke the Maxwell Garnett and Bruggeman formahsms. In fact, the 
bilocally-approximated SPFT gives rise to the Bruggeman homogenization formalism in the 
limit L ^ [12]. 

In the following sections, we consider Voigt wave propagation in a biaxial two-phase 
HCM within the context of the SPFT. A generalized SPFT implementation is developed in 
order to explore the propagation of Voigt waves in any direction. The unit Cartesian vectors 
are denoted as x, y and z. Double underlined quantities are 3x3 dyadics. The wavenumber 
of free space (i.e., vacuum) is /cq. 

2 Homogenization background 

2.1 Component phases 

The propagation of Voigt waves in a two-phase homogenized composite medium (HCM) 
is investigated. Both component material phases are taken as uniaxial dielectric mediums. 
Therefore, they do not individually support Voigt wave propagation [7]. 

Let the components phases — labelled as a and h — be characterized by the relative 
permittivity dyadics 

^=&M ' [e^xx + e^ (yy + zz)] • ^^(99) 
k = CfeXx + efe (yy + zz) 

respectively. The preferred axis of component material phase a is rotated under the action 
of 

R^{(p) = cos<^(xx + yy) +sin9?(xy-yx) +ZZ, (2) 

to lie in the xy plane at an angle ip to the x axis, whereas the preferred axis of component 
phase b is aligned with the x axis, without loss of generality. The superscript indicates the 
transpose operation. 

Let the regions occupied by component phases a and b be denoted by Va and H, respec- 
tively. The component phases are randomly distributed such that all space V — Va U Vf,. 
Spherical microstructural geometries, with characteristic length scales which are small in 
comparison with electromagnetic wavelengths, are assumed for both component phases. The 
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distributional statistics of the component phases are described in terms of moments of the 
characteristic functions 



'1, rey,, 

$,(r)= (e^a,b). (3) 

The volume fraction ff of phase £ is given by the first statistical moment of i.e., ( ^^(r) ) = 
fi . Clearly, fa + fb = 1- The second statistical moment of provides a two-point covariance 
function. We adopt the physically motivated form [13] 

r ($,(r))($,(r')), |r-r'|>L 
($,(r)$,(r'))= , (4) 

[ ($,(r)), |r-r'| <L 

wherein L > is the correlation length. Thus, coherent interactions between a scattering 
centre located at r and another located at r' are accommodated provided that |r — r'| < L. 
However, if |r — r'| > L then the scattering centre at r is presumed to act independently 
of the scattering centre at r'. In implementations of the SPFT, the precise form of the 
covariance function is relatively unimportant to estimate the constitutive parameters of the 
HCM [14]. 



2.2 Homogenized composite medium 

Since the preferred axes of the uniaxial component phases are not generally aligned, the 
HCM is a biaxial dielectric medium. We confine ourselves to the bilocally approximated 
SPFT in order to estimate the HCM relative permittivity dyadic 

ifta = 4a^^ + 4ayy + 4a^^ + 4a{^y + y^)- (5) 

The SPFT is based upon iterative refinements of a comparison medium. The relative per- 
mittivity dyadic of the comparison medium, namely, 

^Br ^ e^^xx + e^^yy + eB^zz + eB^ (xy + yx), (6) 

is provided by the Bruggeman homogenization formalism [12]. 



2.3 Depolarization and polarizability dyadics 

The depolarization dyadic D is central to both the Bruggeman formalism and the SPFT. 
It provides the electromagnetic response of a infinitesimally small spherical exclusion re- 
gion, immersed in a homogeneous background. For the comparison medium with relative 
permittivity dyadic f the corresponding depolarization dyadic is given by [15] 

D = -V r d<P r d(^—6 — ^ qq^ (7) 

= iujAt: Jo ^ Jo rl sin2 ^ + el, cos2 ^ 
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wherein 

^Br = ^Br cos^ + sin^ (f) + 2e^^ sin (f) cos 0, (8) 

and q = sin 9 cos x + sin 9 sin y + cos ^ z is the unit position vector. 

A related construction, much used in homogenization formahsms, is the polarizabihty 
density dyadic x . It is defined here as 



{u-^=Br) ' [L + tujR- (i.-i^J] \ {£ = a,b). (9) 



2.4 The bilocally approximated SPFT 

After accommodating higher-order distributional statistics, the bilocally approximated SPFT 
estimate 

^,a=^Br-^U + ^^a'Ry' ' (1°) 

is derived [12]. The mass operator [16] term 

&. = (2.-xJ -fi-la.-aj (11) 

is specified in terms of the principal value integral 

P = V I d^RGg^(R), (12) 

— Ja<L — 

with R = r — r' and ^^^(R) being the unbounded dyadic Green function of the comparison 
medium. A surface integral representation of E_ is established in the Appendix. Thereby, 
we see that e^^ has a complex dependency upon the correlation length L, with e^^ becoming 
equal to in the limit L — > 0. 



3 Voigt wave propagation 

In order to explore the possibility of Voigt wave propagation — in any direction — it is 
expedient to consider the HCM relative permittivity dyadic (5) in the rotated coordinate 
frame specified by the three Euler angles a, (5 and 7 [17]. On rotation, transforms to 

7) = 1.(7) • ^yi^) ' ^-M) ' Ua ' • • £(7) (13) 

= eiixx + e22yy + e33zz + ei2 (xy + yx) 

+ei3 (xz + zx) + 623 (yz + zy), (14) 

wherein 

RyiP) — cos /9(xx + zz)+sin/?(zx — xz) + yy. (15) 
Two conditions must be satisfied in order for Voigt waves to propagate along the z axis 

[7]: 
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(i) Y{a, P, 7) = 0, and 

(ii) W{a,P,j)j^O. 

In terms of the relative permittivity scalars on the right side of (14) (which are all implicit 
functions of a, (3 and 7), 

y (a, /3, 7) = 613 + 623 - 2623633 [ 2612613 - ( 611 - 622 ) 623 ] + [( ^ll " ^22 f + ] ^33 

+2613 { £23^13 - [ 2612623 + ( 611 - 622 ) 613 ] 633 } (16) 

and 

1^(0;, /5, 7) = 612633 - 613623 • (17) 

4 Numerical results 

The numerical calculations proceed in two stages: Firstly, is estimated using the bilo- 
cally approximated SPFT for a representative example; secondly, the quantities y(a,/9, 7) 
and M^(a,/3, 7) are calculated as functions of the Euler angles. In particular, the (a,/?, 7) 
angular coordinates of the zeros of F, and the corresponding values of W at those (q;,/9, 7) 
coordinates, are sought. The 7 angular coordinate is disregarded since propagation parallel 
to the z axis (of the rotated coordinate system) is independent of rotation about that axis. 

The following constitutive parameters were selected for the component phases a and h 
for all results presented here: 



(18) 



6« = 2.5 + i0.l5, 6« = 6.8 + i0.25 5 
e& = 3.6 + i2.55, 6^-2.2 + il.45 

with the dissipation parameter b e [0,1]. The volume fraction — 0.3 for all calculations. 
4.1 HCM constitutive parameters 

Consider the four relative permittivity scalars, e^^^'^'*, in the unrotated reference frame, i.e., 
6^^(0,0,0) = 6^^^. The 6^^^'^'* values calculated with the dissipation parameter 5 = 1 are 
plotted in Figure 1, as functions of the orientation angle of component phase a and the 
relative correlation length L/cq. At 93 = (and also at </? = tt), the preferred axes of both 
component phases are ahgned. Accordingly, the HCM is uniaxial with 6^„ = 6^„. For </? 7^ 
(and also <^ 7^ tt), the HCM is biaxial. As — > 7r/2, the HCM biaxial structure becomes 
orthorhombic since 6^^ 0. For intermediate values of (/9 G (0,7r/2), the HCM has the 
general non-orthorhombic biaxial form [18]. The correlation length is found to have only a 
marginal influence on e^^'^'^ for 5 — 1. 

The HCM constitutive parameters corresponding to those of Figure 1, but arising from 
nondissipative component phases (i.e., 6 = 0), are presented in Figure 2. The absence of 
dissipation in the component phases has little effect on the real parts of e^^'^'*. However, the 
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imaginary parts of e^^^'^'* are much altered. Since the component phases are nondissipative, 
the imaginary parts of e^^^'^'* are null-valued at zero correlation length. As the correlation 
length increases, the loss due to the effects of coherent scattering becomes greater. Hence, 
the magnitudes of the imaginary parts of e^^^'^' are observed to increase in Figure 2 as Lko 
grows. Furthermore, it is clear from Figure 2 that the rate of increase of these imaginary 
parts is sensitively dependent upon the orientation angle (p of the component phase a. 

4.2 Zeros of Y 

The condition y = can be satisfied at two distinct HCM orientations. These orientations 
are denoted by the angular coordinates (ai. Pi) and (0:2, Pt)- With the normalized correlation 
length fixed at LA;,, = 0.1, the {ai, j3i) and {a2'! P2) angular coordinates arc graphed as 
functions of (/? in Figure 3 for the dissipation parameter values 6 = 0.0,0.1,0.5 and 1. In 
particular, observe that the two distinct solutions of Y{a,P,0) = exist even when the 
component material phases are nondissipative (i.e., 5 — 0). The angular coordinates {ai, (3i) 
and {a2, (32) are clearly sensitive to both (f and 5. 

Values of \W\, corresponding to the angular coordinates {ai, f3i) and {a2, ^2) of Figure 3, 
are plotted against (/? in Figure 4. For (/? > 0, the magnitude \W\ > 0. In particular, the 
inequality > holds for 5 = (which is not clearly illustrated in Figure 4 due to limited 
resolution). Therefore, Voigt waves can propagate along two distinct singular axes in the 
biaxial HCM, as specified by the angular coordinates {ai, and (0:2, /32), even when the 
HCM arises from nondissipative component phases. This conclusion stems solely from the 
incorporation of the correlation length in the SPFT, because the Maxwell Garnett and the 
Bruggeman formalisms would not predict Voigt wave propagation when both component 
phases are nondissipative [2]. 

The two orientations that zero the value of Y, as specified by the angular coordinates 
{ai, (3i) and («2i/32), are plotted against (/? in Figure 5 for the normalized correlation lengths 
Lko = 0,0.05 and 0.1. The dissipation parameter is fixed at 5 = 1. As in Figure 3, the 
two distinct directions described by {ai, Pi) and {a2, P2) are sensitively dependent upon the 
orientation angle ip of the component phase a. Furthermore, the two distinct directions 
persist in the limit Lko 0. The infiuence of Lko upon the angular coordinates (cci. Pi) and 
{0:2, P2) (as illustrated in Figure 5) is relatively minor in comparison with the infiuence of 
the dissipation parameter 6 (as illustrated in Figure 3). 

For the angular coordinates (ai,/3i) and {a2, P2) of Figure 5, the corresponding values 
of \W\ are presented in Figure 6 as functions of </?. Clearly, \W\ > for (/? > when 
Lko — 0,0.05 and 0.1. The magnitude of Lko has only a minor infiuence on \W\. Hence, 
the orientations of the singular axes, along which Voigt waves may propagate in the biaxial 
HCM, are modulated only to a minor degree by the correlation length. 
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5 Conclusions 



The role of the correlation length in facilitating the propagation of Voigt waves in HCMs is 
delineated. Thereby, the importance of taking higher-order distributional statistics in ho- 
mogenization studies into account is further emphasized. Specifically, we have demonstrated 
that 

1. Voigt waves can propagate in HCMs arising from nondissipative component phases, 
provided that a nonzero correlation length is accommodated, according to the SPFT. 

2. The orientations of singular axes in HCMs are sensitively dependent upon (i) the 
degree of dissipation exhibited by the component phases and the (ii) the orientation of 
the preferred axes of the component material phases. By comparison, the correlation 
length plays only a secondary role in determining the singular axis orientations. 
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Appendix 

We establish here a surface integral representation of P (12), amenable to numerical evalu- 
ation. A straightforward specialization of the evaluation of ^ for bianisotropic HCMs [12] 
yields the volume integral 



p fafb 

— 27r2 iuj 



{q/ujfa + (3 



sin qL 



{q/uj)Hc + {q/uj)HB + tA V Q 



L cos qL 



where the scalar quantities tA,tB and tc are given as 



tB = -A^o|eB,(e^, + e^,)cos'^ + 
tc = jjio (e%^ cos^ 6 + sin^ 6^ . 
The dyadic quantities a and (3 are specified as 



p0 z , X y 
^ Br^Br ^ ^Br^Br 



Br 



sm 



a 



= B 



Za- 



tB 


T 


tc 




tA 


T 


tc 





C" 



(19) 

(20) 

(21) 
(22) 

(23) 
(24) 
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with 



T 

= B 



adj (i^ J 



(26) 



T 
=c 



+ 
+ 

+ 

1^0 



'Br 



+ € 



'Br 



-Br 



sm 



sin^^ 



XX 



- cos^ (f)^ sin^ 6 
{el, + el,) cos' e + V%, sin' 



yy 

z z 



Br 



sm cos ( 



'Br 



sin'^ 



(xy + yx) }, 



sin' 9 cos' X X + sin' 9 sin' y y 



+ cos' ^ z z + sin' 9 cos sin (xy + yx) 



Let 



2ic 



(26) 
(27) 

(28) 



In the long-wavelength regime, i.e., | Ly/p± | <^ 1 [12], the application of residue calculus to 
(19) delivers the surface integral 



P = 



Ani Jo 



d9 



sin 9 



3 - At Ate 



X 



2 (P+ 



a + iL(VpT-V^)^ . (29) 



Standard numerical techniques may be straightforwardly applied to evaluate (29). 
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Figure 1: The real (left) and imaginary (right) parts of the components of £^^(0, 0, 0) = e^^ 
plotted against the relative correlation length Lkg and orientation angle ip (in degrees) of 
component phase a. The dissipation parame^^ 6 = 1. 
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Figure 3: The a and /3 angular coordinates (in degrees) of the zeros of Y{a,(3,0) plotted 
against the orientation angle ip (in degrees) of component phase a. The relative correlation 
length Lko = 0.1. 
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Figure 4: The values of |Vr(Q;i,2, /3i,2, 0)| corresponding to the (ai,/5i) and {a2, P2) angular 
coordinates of Figure 3, as functions of the orientation angle (p (in degrees) of the component 
phase a. 
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Figure 6: The values of |VF(q;i,2, /3i,2, 0)| corresponding to the (q;i,/3i) and (q;2,/?2) angular 
coordinates of Figure 5, as functions of the orientation angle </? (in degrees) of the component 
phase a. 
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